module interp

    implicit none

contains

    subroutine interp_linear(nd, n, x, v, x0, v0)
        integer, intent(in) :: nd, n
        real, intent(in) :: x(nd,n)
        real, intent(in) :: v(n)
        real, intent(in) :: x0(nd)
        real, intent(out) :: v0

        real xx, yy, a, b, c, d

        if(nd == 1 .and. n == 2) then
            xx = (x0(1)-x(1,1))/(x(1,2)-x(1,1))
            a = v(1)
            b = v(2)-v(1)
            v0 = a+b*xx
        else if(nd == 2) then
            if(n == 2) then
                xx = sqrt((x0(1) -x(1,1))**2+(x0(2) -x(2,1))**2)/ &
                     sqrt((x(1,1)-x(1,2))**2+(x(2,1)-x(2,2))**2)
                a = v(1)
                b = v(2)-v(1)
                v0 = a+b*xx
            else if(n == 4) then
                xx = (x0(1)-x(1,1))/(x(1,2)-x(1,1))
                yy = (x0(2)-x(2,1))/(x(2,3)-x(2,1))
                a = v(1)
                b = v(2)-v(1)
                c = v(3)-v(1)
                d = v(1)-v(2)-v(3)+v(4)
                v0 = a+b*xx+c**yy+d*xx*yy
            end if
        end if

        return
    end subroutine interp_linear
    
    subroutine interp_linear_derv_x(nd, n, x, v, x0, dv0)
        integer, intent(in) :: nd, n
        real, intent(in) :: x(nd,n)
        real, intent(in) :: v(n)
        real, intent(in) :: x0(nd)
        real, intent(out) :: dv0
        
        if(nd == 1 .and. n == 2) then
            dv0 = (v(1)-v(2))/(x(1,1)-x(1,2))
        else if(nd == 2) then
            if(n == 2) then
                ! ---> under construction
            else if(n == 4) then
                ! ---> under construction
            end if
        end if
        
        
        return
    end subroutine interp_linear_derv_x
    
    subroutine interp_linear_derv_y(nd, n, x, v, x0, dv0)
        integer, intent(in) :: nd, n
        real, intent(in) :: x(nd,n)
        real, intent(in) :: v(n)
        real, intent(in) :: x0(nd)
        real, intent(out) :: dv0
        
        if(nd == 2) then
            if(n == 2) then
                ! ---> under construction
            else if(n == 4) then
                ! ---> under construction
            end if
        end if
        
        return
    end subroutine interp_linear_derv_y
    
    subroutine interp_linear_derv_z(nd, n, x, v, x0, dv0)
        integer, intent(in) :: nd, n
        real, intent(in) :: x(nd,n)
        real, intent(in) :: v(n)
        real, intent(in) :: x0(nd)
        real, intent(out) :: dv0
        
        if(nd == 3) then
            ! ---> under construction
        end if
        
        return
    end subroutine interp_linear_derv_z
    
end module interp
